### PPP Further Exploration 

rm(list = ls())

### libraries and packages
pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", "fracdiff",
         "fractal",
         "lme4", "ArfimaMLM","gam",
         "forecast",
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)
# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

# creating a dummy for more than 150 AQI
p2$aqi_high <- 0
p2$aqi_high[p2$aqi>150] <- 1

###
day_pollute.h <- lm(day_pollute ~ aqi_high  + national+ parade, data=p2)
vc.day_pollute.h <-cluster.vcov(day_pollute.h, p2$day) 
se.vc.day_pollute.h <- coeftest(day_pollute.h, vc.day_pollute.h)

day_pollute.hc <- lm(day_pollute ~ aqi_high  + national + parade + hukou + insider + soe 
+ kids + educ_sd + married + gender + ccp + age_sd + income_sd, data=p2)
vc.day_pollute.hc <-cluster.vcov(day_pollute.hc, p2$day) 
se.vc.day_pollute.hc <- coeftest(day_pollute.hc, vc.day_pollute.h)

## small table for this
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables")

fitlist.high <- list(day_pollute.h,day_pollute.hc)
se.high <- list(se.vc.day_pollute.h[,2],se.vc.day_pollute.hc[,2])
cov.labs <- c("Unhealthy Pollution (Dummy)", "National Holiday", "Parade Period","Hukou Status", "Regime Insider", "SOE Employee", "Children", "Education", "Married", "Gender", "CCP Member", "Age", "Income")
omit.stats <- c("adj.rsq","ser","f")
## TABLE A3
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/high.perception.tex")
stargazer(title = "Unhealthy Pollution Levels and Perception Of Severity",
          se = se.high,
          header = F,
          fitlist.high,
          covariate.labels = cov.labs,
          omit.stat = omit.stats,
          no.space = T,
          suppress.errors = F,float = F,
          omit.table.layout = "n",
          dep.var.labels =c("Perception That Daily Pollution is Not a Serious Problem")
)
sink()

